Ordinal Logistic Regression
Advanced Categorical Data Analysis
This document serves both learning and practical purposes. It is designed for educational use, aiming to enhance statistical analysis skills and provide clear, organized notes for future reference.
1 Introduction
Ordinal logistic regression is a statistical analysis used for predicting an ordinal outcome based on one or more predictor variables. An ordinal outcome has categories with a natural order but unknown distances between categories, such as ratings (e.g., poor, fair, good, excellent).
The dependent variable is ordinal, meaning it has a meaningful order. Examples include survey responses (e.g., strongly disagree to strongly agree) or educational attainment (e.g., high school, bachelor’s, master’s, doctorate).
Ordinal logistic regression includes several models that are used to analyze ordinal data. The most common models are:
- Adjacent-Category Logit Model
- Continuation-Ratio Logit Model
- Proportional Odds Model (Cumulative Logit Model)
1.1 Adjacent-Category Logit Model
This model compares each response category with the next adjacent category. The model can be written as:
\[ \log \left( \frac{P(Y = j)}{P(Y = j+1)} \right) = \alpha_j + \beta X \]
- \(Y\) is the ordinal response variable.
- \(j\) is the category.
- \(\alpha_j\) are the intercepts for the \(j\)-th category.
- \(\beta\) is the coefficient for the predictor variable \(X\).
The model assumes that the log-odds do not depend on the specific category being compared. The model helps understand how the odds change from one category to the next adjacent category based on the predictor variables.
1.2 Continuation-Ratio Logit Model
This model is used to compare each response category with all lower categories combined. The model can be written as:
\[ \log \left( \frac{P(Y = j)}{P(Y < j)} \right) = \alpha_j + \beta X \]
- \(Y\) is the ordinal response variable.
- \(j\) is the category.
- \(\alpha_j\) are the intercepts for the \(j\)-th category.
- \(\beta\) is the coefficient for the predictor variable \(X\).
Different constant terms and slopes are allowed, making it more flexible. The model compares the probability of being in a particular category to the probability of being in all lower categories combined, providing detailed insights into transitions between levels.
1.3 Proportional Odds Model (Cumulative Logit Model)
The Proportional Odds Model is used to model the cumulative probabilities of the ordinal outcome variable. The model compares the cumulative probability of the response being in a category less than or equal to a specific value versus being in a higher category:
\[ \log \left( \frac{P(Y \leq j)}{P(Y > j)} \right) = \alpha_j + \beta X \]
- \(Y\) is the ordinal response variable.
- \(j\) is the category.
- \(\alpha_j\) are the intercepts for the \(j\)-th category.
- \(\beta\) is the coefficient for the predictor variable \(X\).
The relationship between predictors and the log-odds of being in a lower versus higher category is constant across all categories. The coefficients \(\beta\) indicate the effect of the predictors on the log-odds of the outcome being at or below a certain category versus above it.
2 Setting Up the Environment
Load the necessary libraries to ensure the R environment is equipped with the tools and functions required for efficient and effective analysis.
3 Practical 1
The wine dataset used in this exercise is used to illustrate how to apply cumulative link models (Proportional Odds Models) for analyzing ordinal data. This dataset contains 72 observations and 6 variables.
- response: A numeric variable representing the response value.
- rating: An ordinal factor variable with levels indicating the bitterness rating of the wine (from 1 to 5, with 1 being the least bitter and 5 being the most bitter).
- temp: A factor variable indicating the temperature at which the wine was served (e.g., cold, warm).
- contact: A factor variable indicating whether the wine had contact with air (e.g., yes, no).
- bottle: A numeric variable representing different bottle codes.
- judge: A factor variable representing different judges who rated the wine.
3.1 Data Reading and Exploration
First, load the dataset and convert any labeled variables to factors. Then, perform exploratory data analysis to understand the data distribution and relationships between variables.
# Reading the data
data.o1 <- read_dta("Datasets/wine.dta")
# Briefly examine the data
glimpse(data.o1)Rows: 72
Columns: 6
$ response <dbl> 36, 48, 47, 67, 77, 60, 83, 90, 17, 22, 14, 50, 30, 51, 90, 7…
$ rating <dbl+lbl> 2, 3, 3, 4, 4, 4, 5, 5, 1, 2, 1, 3, 2, 3, 5, 4, 2, 3, 3, …
$ temp <dbl+lbl> 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, 1, 2, 2, 2, 2, 1, 1, 1, …
$ contact <dbl+lbl> 1, 1, 2, 2, 1, 1, 2, 2, 1, 1, 2, 2, 1, 1, 2, 2, 1, 1, 2, …
$ bottle <dbl+lbl> 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, …
$ judge <dbl+lbl> 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, …
# Summarize the data
summary(data.o1) response rating temp contact bottle
Min. :12.00 Min. :1.000 Min. :1.0 Min. :1.0 Min. :1.00
1st Qu.:32.00 1st Qu.:2.000 1st Qu.:1.0 1st Qu.:1.0 1st Qu.:2.75
Median :46.00 Median :3.000 Median :1.5 Median :1.5 Median :4.50
Mean :47.22 Mean :2.917 Mean :1.5 Mean :1.5 Mean :4.50
3rd Qu.:60.00 3rd Qu.:4.000 3rd Qu.:2.0 3rd Qu.:2.0 3rd Qu.:6.25
Max. :90.00 Max. :5.000 Max. :2.0 Max. :2.0 Max. :8.00
judge
Min. :1
1st Qu.:3
Median :5
Mean :5
3rd Qu.:7
Max. :9
# Convert labeled variables to factor variables
data.o2 <- data.o1
data.o2 <- data.o2 %>% mutate(across(where(is.labelled), ~ as_factor(.x)))
# Check the data structure after conversion
glimpse(data.o2)Rows: 72
Columns: 6
$ response <dbl> 36, 48, 47, 67, 77, 60, 83, 90, 17, 22, 14, 50, 30, 51, 90, 7…
$ rating <fct> 2, 3, 3, 4, 4, 4, 5, 5, 1, 2, 1, 3, 2, 3, 5, 4, 2, 3, 3, 2, 5…
$ temp <fct> cold, cold, cold, cold, warm, warm, warm, warm, cold, cold, c…
$ contact <fct> no, no, yes, yes, no, no, yes, yes, no, no, yes, yes, no, no,…
$ bottle <fct> 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5, 6, 7, 8, 1, 2, 3, 4, 5…
$ judge <fct> 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 3…
# Summarize the characteristics of the dataset by `rating`
data.o2 %>%
select(response, temp, rating, contact) %>%
tbl_summary(by = rating)| Characteristic | 1, N = 51 | 2, N = 221 | 3, N = 261 | 4, N = 121 | 5, N = 71 |
|---|---|---|---|---|---|
| response | 14 (13, 17) | 31 (27, 35) | 48 (45, 51) | 67 (62, 72) | 84 (82, 89) |
| temp | |||||
| cold | 5 (100%) | 16 (73%) | 13 (50%) | 2 (17%) | 0 (0%) |
| warm | 0 (0%) | 6 (27%) | 13 (50%) | 10 (83%) | 7 (100%) |
| contact | 1 (20%) | 8 (36%) | 13 (50%) | 9 (75%) | 5 (71%) |
| 1 Median (IQR); n (%) | |||||
3.2 Estimation
The objective is to estimate the cumulative link model (proportional odds model) using the clm function from the ordinal package. This involves estimating regression coefficients, standard errors, and p-values. The cumulative link model (logit link) is specified as:
\[ P(Y_i \leq j) = \theta_j - \beta_1(\text{temp}) - \beta_2(\text{contact}) \]
\(\theta_j\) are the threshold parameters, and \(\beta_1\) and \(\beta_2\) are the regression coefficients for temp and contact, respectively.
# Fit the cumulative link model
mod.o1 <- clm(rating ~ temp + contact, data = data.o2)
summary(mod.o1)formula: rating ~ temp + contact
data: data.o2
link threshold nobs logLik AIC niter max.grad cond.H
logit flexible 72 -86.49 184.98 6(0) 4.02e-12 2.7e+01
Coefficients:
Estimate Std. Error z value Pr(>|z|)
tempwarm 2.5031 0.5287 4.735 2.19e-06 ***
contactyes 1.5278 0.4766 3.205 0.00135 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Threshold coefficients:
Estimate Std. Error z value
1|2 -1.3444 0.5171 -2.600
2|3 1.2508 0.4379 2.857
3|4 3.4669 0.5978 5.800
4|5 5.0064 0.7309 6.850
Summary of Output
-
Coefficients:
tempwarm: The estimated coefficient for
tempwarmis 2.5031, meaning that having a warm temperature increases the log odds of being in a higher rating category by 2.5031, compared to the baseline category (cold temperature). The effect is statistically significant at the 0.001 level (p-value = 2.19e-06).contactyes: The estimated coefficient for
contactyesis 1.5278, meaning that having contact increases the log odds of being in a higher rating category by 1.5278, compared to the baseline category (no contact). The effect is statistically significant at the 0.01 level (p-value = 0.00135).
-
Threshold Coefficients
1|2: The threshold coefficient of -1.3444 separates the first category from the second. This value indicates the log odds of being in category 1 or lower versus category 2 or higher.
2|3: The threshold coefficient of 1.2508 separates the second category from the third. This value indicates the log odds of being in category 2 or lower versus category 3 or higher.
3|4: The threshold coefficient of 3.4669 separates the third category from the fourth. This value indicates the log odds of being in category 3 or lower versus category 4 or higher.
4|5: The threshold coefficient of 5.0064 separates the fourth category from the fifth. This value indicates the log odds of being in category 4 or lower versus category 5.
-
Interpretation
tempwarm: The estimate is 2.5031 with a p-value of 2.19e-06, indicating a significant positive association between warmer temperature and higher bitterness ratings.
contactyes: The estimate is 1.5278 with a p-value of 0.00135, indicating a significant positive association between contact with air and higher bitterness ratings.
# Tidy summary of the model
tidy(mod.o1)# A tibble: 6 × 6
term estimate std.error statistic p.value coef.type
<chr> <dbl> <dbl> <dbl> <dbl> <chr>
1 1|2 -1.34 0.517 -2.60 9.33e- 3 intercept
2 2|3 1.25 0.438 2.86 4.28e- 3 intercept
3 3|4 3.47 0.598 5.80 6.64e- 9 intercept
4 4|5 5.01 0.731 6.85 7.41e-12 intercept
5 tempwarm 2.50 0.529 4.73 2.19e- 6 location
6 contactyes 1.53 0.477 3.21 1.35e- 3 location
Interpretation
1|2: The threshold separating category 1 from category 2 has a log-odds estimate of -1.34, which is statistically significant (p-value < 0.01).
2|3: The threshold separating category 2 from category 3 has a log-odds estimate of 1.25, which is statistically significant (p-value < 0.01).
3|4: The threshold separating category 3 from category 4 has a log-odds estimate of 3.47, which is statistically significant (p-value < 0.001).
4|5: The threshold separating category 4 from category 5 has a log-odds estimate of 5.01, which is statistically significant (p-value < 0.001).
tempwarm: The log-odds of the response variable being in a higher category increase by 2.50 when the temperature is warm compared to cold. This effect is statistically significant (p-value < 0.001).
contactyes: The log-odds of the response variable being in a higher category increase by 1.53 when there is contact compared to no contact. This effect is statistically significant (p-value < 0.01).
# Tidy summary with exponentiated coefficients and confidence intervals
tidy(mod.o1, exponentiate = TRUE, conf.int = TRUE)# A tibble: 6 × 8
term estimate std.error statistic p.value conf.low conf.high coef.type
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr>
1 1|2 0.261 0.517 -2.60 9.33e- 3 NA NA intercept
2 2|3 3.49 0.438 2.86 4.28e- 3 NA NA intercept
3 3|4 32.0 0.598 5.80 6.64e- 9 NA NA intercept
4 4|5 149. 0.731 6.85 7.41e-12 NA NA intercept
5 tempwarm 12.2 0.529 4.73 2.19e- 6 4.53 36.4 location
6 contactyes 4.61 0.477 3.21 1.35e- 3 1.85 12.1 location
Interpretation
1|2: The odds of being in category 1 or lower versus category 2 or higher are 0.261 times, indicating a decrease in odds, which is statistically significant.
2|3: The odds of being in category 2 or lower versus category 3 or higher are 3.49 times, indicating an increase in odds, which is statistically significant.
3|4: The odds of being in category 3 or lower versus category 4 or higher are 32 times, indicating a substantial increase in odds, which is statistically significant.
4|5: The odds of being in category 4 or lower versus category 5 are 149 times, indicating a substantial increase in odds, which is statistically significant.
tempwarm: The odds of the response variable being in a higher category are 12.2 times higher when the temperature is warm compared to cold. This effect is statistically significant, and the confidence interval (4.53 to 36.4) indicates the range within which the true odds ratio is likely to fall with 95% confidence.
contactyes: The odds of the response variable being in a higher category are 4.61 times higher when there is contact compared to no contact. This effect is statistically significant, and the confidence interval (1.85 to 12.1) indicates the range within which the true odds ratio is likely to fall with 95% confidence.
3.3 Inferences
The inferences from the fitted cumulative link model are performed by:
- Checking Wald-based p-values
- Calculating confidence intervals
- Comparing models
3.3.1 Wald-Based P-Values
Wald-based p-values test whether each parameter (regression coefficient) is significantly different from zero. These were already included in the summary and tidy outputs, showing which predictors are significant.
3.3.2 Calculating Confidence Intervals
Confidence intervals provide a range of values within which the true parameter value is likely to fall. The confidence intervals were displayed in the previous tidy output with conf.int = TRUE.
3.3.3 Model Comparison using Likelihood Ratio Test
Let’s compare two nested models using the likelihood ratio test:
-
Model
mod.o1b: A simpler model with onlytempas the predictor. -
Model
mod.o1: A more complex model including bothtempandcontactas predictors.
# Fit the simpler model
mod.o1b <- clm(rating ~ temp, data = data.o2)
# Perform the likelihood ratio test to compare models
anova(mod.o1b, mod.o1, test = "Chisq")'test' argument ignored in anova.clm
Likelihood ratio tests of cumulative link models:
formula: link: threshold:
mod.o1b rating ~ temp logit flexible
mod.o1 rating ~ temp + contact logit flexible
no.par AIC logLik LR.stat df Pr(>Chisq)
mod.o1b 5 194.03 -92.013
mod.o1 6 184.98 -86.492 11.043 1 0.0008902 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Summary of Output
The p-value of 0.0008902 is highly significant, indicating that adding
contactto the model significantly improves the fit compared to the simpler model with onlytemp.Therefore, model
o1(which includes bothtempandcontact) is preferred over modelo1b.
3.4 Checking Proportional Odds Assumption
The proportional odds assumption is crucial for the cumulative link model. It implies that the relationship between each pair of outcome groups is the same. There are two methods to check the proportional odds assumption.
3.4.1 Method 1: Using nominal_test Function
The ordinal::nominal_test() function tests the ordinality of each predictor. If the p-value is greater than 0.05, it fails to reject the null hypothesis that the proportional odds assumption holds for that predictor.
nominal_test(mod.o1)Tests of nominal effects
formula: rating ~ temp + contact
Df logLik AIC LRT Pr(>Chi)
<none> -86.492 184.98
temp 3 -84.904 187.81 3.1750 0.3654
contact 3 -86.209 190.42 0.5667 0.9040
Interpretation:
For both
tempandcontact, the p-values are greater than 0.05 (0.3654 and 0.9040, respectively).This indicates that the assumption of proportional odds cannot be rejected for these variables, suggesting that the proportional odds assumption holds.
3.4.2 Method 2: Likelihood Ratio Test
This method involves comparing two models:
-
Model
mod.o1- This model includes the predictors
tempandcontactand assumes that the proportional odds assumption holds for all predictors.
- This model includes the predictors
-
Model
mod.o1.nominal:This model treats the predictor
contactas having nominal effects, meaning it does not assume proportional odds for this predictor.Instead, it allows the effect of
contactto vary across different thresholds of the ordinal response variable.
# Fit the nominal model
mod.o1.nominal <- clm(rating ~ temp, nominal = ~contact, data = data.o2)nominal = ~contact specifies that contact should be treated as having non-proportional effects.
# Compare the two models using `anova` function
anova(mod.o1, mod.o1.nominal)Likelihood ratio tests of cumulative link models:
formula: nominal: link: threshold:
mod.o1 rating ~ temp + contact ~1 logit flexible
mod.o1.nominal rating ~ temp ~contact logit flexible
no.par AIC logLik LR.stat df Pr(>Chisq)
mod.o1 6 184.98 -86.492
mod.o1.nominal 9 190.42 -86.209 0.5667 3 0.904
Interpretation:
The p-value is not significant at the 5% level, indicating that the more complex model (
o1.nominal) does not provide a significantly better fit than the simpler model (o1).The proportional odds assumption holds for the predictor
contact. Therefore, the simpler model (o1), which assumes proportional odds forcontact, is adequate.
3.5 Prediction
3.5.1 Method 1: Predicting Probabilities for Existing Data
This method focuses on predicting the probabilities (fitted probabilities) that each observation in the existing dataset falls into each response category.
The outcome variable (rating) is removed from the dataset because it is the response variable we want to predict.
# Predict the probabilities for each category
prob.mod.o1 <- predict(mod.o1, newdata = new.data1)
head(prob.mod.o1$fit) 1 2 3 4 5
1 0.20679013 0.5706497 0.1922909 0.02361882 0.00665041
2 0.20679013 0.5706497 0.1922909 0.02361882 0.00665041
3 0.05354601 0.3776461 0.4430599 0.09582084 0.02992711
4 0.05354601 0.3776461 0.4430599 0.09582084 0.02992711
5 0.02088771 0.2014157 0.5015755 0.20049402 0.07562701
6 0.02088771 0.2014157 0.5015755 0.20049402 0.07562701
tail(prob.mod.o1$fit) 1 2 3 4 5
67 0.053546010 0.37764614 0.4430599 0.09582084 0.02992711
68 0.053546010 0.37764614 0.4430599 0.09582084 0.02992711
69 0.020887709 0.20141572 0.5015755 0.20049402 0.07562701
70 0.020887709 0.20141572 0.5015755 0.20049402 0.07562701
71 0.004608274 0.05380128 0.3042099 0.36359581 0.27378469
72 0.004608274 0.05380128 0.3042099 0.36359581 0.27378469
The rows correspond to different observations in the dataset.
The columns represent the predicted probabilities for each response category (1 to 5).
$fit
[1] 2 2 3 3 3 3 4 4 2 2 3 3 3 3 4 4 2 2 3 3 3 3 4 4 2 2 3 3 3 3 4 4 2 2 3 3 3 3
[39] 4 4 2 2 3 3 3 3 4 4 2 2 3 3 3 3 4 4 2 2 3 3 3 3 4 4 2 2 3 3 3 3 4 4
Levels: 1 2 3 4 5
The
predictfunction withtype = 'class'is used to predict the most likely category for each observation.This output shows the predicted most likely category for each observation, represented by the numbers 1 to 5.
3.5.2 Method 2: Predicting Probabilities for New Data
This method involves predicting the probabilities that new data falls into each response category. This is useful for understanding model behavior under different conditions or for making predictions for hypothetical scenarios. Remember the categories of the outcome are:
levels(data.o2$rating)[1] "1" "2" "3" "4" "5"
The linear predictor and predicted probabilites for these categories:
- temp = cold, contact = no
- temp = warm, contact = no
- temp = cold, contact = yes
- temp = warm, contact = yes
# Create a new dataset that includes all combinations of the levels of the predictors
new.data2 <- expand.grid(temp = levels(data.o2$temp), contact = levels(data.o2$contact))
new.data2 temp contact
1 cold no
2 warm no
3 cold yes
4 warm yes
- This creates a new data frame with all combinations of
tempandcontact.
# Get linear predictors for the new data
lp.new.data2 <- predict(mod.o1, newdata = new.data2, type = "linear.predictor")
# Get probabilities for the new data
prob.new.data2 <- predict(mod.o1, newdata = new.data2, type = "prob")
# View the results
lp.new.data2$eta1
1 2 3 4 5
1 -1.344383 1.2508088 3.4668869 5.0064042 100000.00
2 -3.847485 -1.2522932 0.9637849 2.5033022 99997.50
3 -2.872181 -0.2769889 1.9390893 3.4786065 99998.47
4 -5.375283 -2.7800909 -0.5640127 0.9755045 99995.97
$eta2
1 2 3 4 5
1 -100000.0 -1.344383 1.2508088 3.4668869 5.0064042
2 -100002.5 -3.847485 -1.2522932 0.9637849 2.5033022
3 -100001.5 -2.872181 -0.2769889 1.9390893 3.4786065
4 -100004.0 -5.375283 -2.7800909 -0.5640127 0.9755045
prob.new.data2$fit
1 2 3 4 5
1 0.206790132 0.57064970 0.1922909 0.02361882 0.00665041
2 0.020887709 0.20141572 0.5015755 0.20049402 0.07562701
3 0.053546010 0.37764614 0.4430599 0.09582084 0.02992711
4 0.004608274 0.05380128 0.3042099 0.36359581 0.27378469
$eta1 and $eta2: Linear predictors for each category.
$fit: Predicted probabilities for each category for the new data.
# Combine the new data with the predicted probabilities for a comprehensive view
cbind(new.data2, prob.new.data2) temp contact fit.1 fit.2 fit.3 fit.4 fit.5
1 cold no 0.206790132 0.57064970 0.1922909 0.02361882 0.00665041
2 warm no 0.020887709 0.20141572 0.5015755 0.20049402 0.07562701
3 cold yes 0.053546010 0.37764614 0.4430599 0.09582084 0.02992711
4 warm yes 0.004608274 0.05380128 0.3042099 0.36359581 0.27378469
3.5.3 Method 3: Using polr Function from MASS Package
This method involves using the polr function from the MASS package to fit a proportional odds model and then predict probabilities for each category of the ordinal response variable.
Re-fitting to get Hessian
Call:
MASS::polr(formula = rating ~ temp + contact, data = data.o2)
Coefficients:
Value Std. Error t value
tempwarm 2.503 0.5287 4.735
contactyes 1.528 0.4766 3.205
Intercepts:
Value Std. Error t value
1|2 -1.3444 0.5171 -2.5998
2|3 1.2508 0.4379 2.8565
3|4 3.4669 0.5978 5.7998
4|5 5.0064 0.7309 6.8496
Residual Deviance: 172.9838
AIC: 184.9838
The model is fitted using the
polrfunction from theMASSpackage. The output provides the estimated coefficients for the predictors and the thresholds (intercepts) for the ordinal response variable.tempwarm: The coefficient for warm temperature, indicating a significant positive association with higher ratings.contactyes: The coefficient for contact with air, indicating a significant positive association with higher ratings.
# Predict the probabilities for each category
prob_polr <- predict(mod_polr, type = "probs")
head(prob_polr) 1 2 3 4 5
1 0.2067917 0.5706466 0.1922920 0.02361917 0.006650532
2 0.2067917 0.5706466 0.1922920 0.02361917 0.006650532
3 0.0535471 0.3776458 0.4430586 0.09582112 0.029927296
4 0.0535471 0.3776458 0.4430586 0.09582112 0.029927296
5 0.0208885 0.2014185 0.5015746 0.20049218 0.075626257
6 0.0208885 0.2014185 0.5015746 0.20049218 0.075626257
tail(prob_polr) 1 2 3 4 5
67 0.053547101 0.37764585 0.4430586 0.09582112 0.02992730
68 0.053547101 0.37764585 0.4430586 0.09582112 0.02992730
69 0.020888502 0.20141847 0.5015746 0.20049218 0.07562626
70 0.020888502 0.20141847 0.5015746 0.20049218 0.07562626
71 0.004608507 0.05380284 0.3042139 0.36359456 0.27378016
72 0.004608507 0.05380284 0.3042139 0.36359456 0.27378016
The
predictfunction withtype = 'probs'is used to obtain the predicted probabilities for each response category for each observation in the dataset.The rows correspond to different observations in the dataset.
The columns represent the predicted probabilities for each response category (1 to 5).
3.5.4 Manual Calculation of Predictions
Manually calculate the probabilities for an ordinal outcome using the coefficients from a cumulative link model mod.o1.
3.5.4.1 Model Summary of mod.o1
Review the model summary to get the coefficients for the predictors and the threshold values needed for calculations.
# Model summary
summary(mod.o1)formula: rating ~ temp + contact
data: data.o2
link threshold nobs logLik AIC niter max.grad cond.H
logit flexible 72 -86.49 184.98 6(0) 4.02e-12 2.7e+01
Coefficients:
Estimate Std. Error z value Pr(>|z|)
tempwarm 2.5031 0.5287 4.735 2.19e-06 ***
contactyes 1.5278 0.4766 3.205 0.00135 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Threshold coefficients:
Estimate Std. Error z value
1|2 -1.3444 0.5171 -2.600
2|3 1.2508 0.4379 2.857
3|4 3.4669 0.5978 5.800
4|5 5.0064 0.7309 6.850
3.5.4.2 New Data for Prediction
The new data points for which predictions are needed are defined using the previously generated new data from Prediction: Method 2.
new.data2 temp contact
1 cold no
2 warm no
3 cold yes
4 warm yes
3.5.4.3 Obtain Linear Predictors
Use the predict function with type = "linear.predictor" to get the linear predictors for the new data points, represented by eta1 and eta2 in the output. Refer Prediction: Method 2.
lp.new.data2$eta1
1 2 3 4 5
1 -1.344383 1.2508088 3.4668869 5.0064042 100000.00
2 -3.847485 -1.2522932 0.9637849 2.5033022 99997.50
3 -2.872181 -0.2769889 1.9390893 3.4786065 99998.47
4 -5.375283 -2.7800909 -0.5640127 0.9755045 99995.97
$eta2
1 2 3 4 5
1 -100000.0 -1.344383 1.2508088 3.4668869 5.0064042
2 -100002.5 -3.847485 -1.2522932 0.9637849 2.5033022
3 -100001.5 -2.872181 -0.2769889 1.9390893 3.4786065
4 -100004.0 -5.375283 -2.7800909 -0.5640127 0.9755045
3.5.4.4 Extract Coefficients
Extract the coefficients from the model that include the estimates for the thresholds and the predictor variables:
- Thresholds: 1|2, 2|3, 3|4, 4|5
- Predictors: tempwarm, contactyes
coef.mod.o1 <- coef(mod.o1)
coef.mod.o1 1|2 2|3 3|4 4|5 tempwarm contactyes
-1.344383 1.250809 3.466887 5.006404 2.503102 1.527798
3.5.4.5 Calculate Linear Predictors for Each Observation
For each new observation, calculate the linear predictors (bx) using the coefficients and the new data values. Using fourth observation as an example:
lp.mod.o1.bx4 <- coef.mod.o1["tempwarm"] * 1 + coef.mod.o1["contactyes"] * 1
lp.mod.o1.bx4tempwarm
4.0309
This calculates the linear predictor for the fourth observation, where temp is warm (coded as 1) and contact is yes (coded as 1).
3.5.4.6 Calculate Logits
Using the linear predictors (bx) and the thresholds, calculate the logits using the formula:
\[ \text{logit}(P(Y \leq j)) = \beta_{0j} - (\beta_1 X_1 + \beta_2 X_2 + \cdots + \beta_p X_p) \]
- \(\beta_{0j}\): Threshold coefficient for category \(j\).
- \(\beta_1, \beta_2, \ldots, \beta_p\): Coefficients for the predictor variables.
- \(X_1, X_2, \ldots, X_p\): Values of the predictor variables.
- Logit Calculation: The linear predictor (\(\beta_1 X_1 + \beta_2 X_2 + \cdots + \beta_p X_p\)) is subtracted from the threshold coefficient to get the logit.
Given the threshold coefficients from the model summary:
- \(\beta_{01|2} = -1.3444\)
- \(\beta_{02|3} = 1.2508\)
- \(\beta_{03|4} = 3.4669\)
- \(\beta_{04|5} = 5.0064\)
logit1 <- coef.mod.o1[1] - lp.mod.o1.bx4
logit2 <- coef.mod.o1[2] - lp.mod.o1.bx4
logit3 <- coef.mod.o1[3] - lp.mod.o1.bx4
logit4 <- coef.mod.o1[4] - lp.mod.o1.bx4
logit1 1|2
-5.375283
logit2 2|3
-2.780091
logit3 3|4
-0.5640127
logit4 4|5
0.9755045
3.5.4.7 Convert Logits to Probabilities
The probabilities represent the cumulative probabilities for each threshold. The logistic function is used to convert logits (log-odds) into probabilities. The formula for the logistic function:
\[ \text{logit}^{-1}(x) = \frac{1}{1 + e^{-x}} \]
\(x\) is the logit value.
3.5.4.8 Calculate Category Probabilities
To find the probability of each specific category, calculate the differences between the cumulative probabilities. This will give the probability of the response falling into each category.
pMat <- cbind(
p1 = pLeq1,
p2 = pLeq2 - pLeq1,
p3 = pLeq3 - pLeq2,
p4 = pLeq4 - pLeq3,
p5 = 1 - pLeq4
)
pMat p1 p2 p3 p4 p5
1|2 0.004608274 0.05380128 0.3042099 0.3635958 0.2737847
3.5.4.9 Confirm Predictions with the Model
Finally, confirm the manual calculations by comparing them with the predictions from the predict function using type = "prob". Refer Method 2.
prob.new.data2$fit
1 2 3 4 5
1 0.206790132 0.57064970 0.1922909 0.02361882 0.00665041
2 0.020887709 0.20141572 0.5015755 0.20049402 0.07562701
3 0.053546010 0.37764614 0.4430599 0.09582084 0.02992711
4 0.004608274 0.05380128 0.3042099 0.36359581 0.27378469
4 Practical 2
The lowbwt.dta dataset used in this exercise is sourced from the Applied Logistic Regression book and involves data that are related to logistic regression models for ordinal outcomes.
Variables:
-
id: Identification number for each observation. -
lbw: Indicator variable for low birth weight. -
age: Age of the mother in years. -
lwt: Weight of the mother at the last menstrual period in pounds. -
race: Race of the mother. -
smoke: Smoking status during pregnancy. -
ptl: Number of previous premature labors. -
hyper: Presence of hypertension. -
urirr: Presence of uterine irritability. -
pvft: Physician visits during the first trimester. -
weight: Birth weight in grams. -
agecat: Age category of the mother. -
wcat: Weight category of the mother. -
anyptl: Indicator variable for any previous premature labor. -
newpvft: New variable for physician visits during the first trimester.
4.1 Data Reading and Exploration
id lbw age lwt
Min. : 4.0 Min. :0.0000 Min. :14.00 Min. : 80.0
1st Qu.: 68.0 1st Qu.:0.0000 1st Qu.:19.00 1st Qu.:110.0
Median :123.0 Median :0.0000 Median :23.00 Median :121.0
Mean :121.1 Mean :0.3122 Mean :23.24 Mean :129.8
3rd Qu.:176.0 3rd Qu.:1.0000 3rd Qu.:26.00 3rd Qu.:140.0
Max. :226.0 Max. :1.0000 Max. :45.00 Max. :250.0
race smoke ptl hyper
Min. :1.000 Min. :0.0000 Min. :0.0000 Min. :0.00000
1st Qu.:1.000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.00000
Median :1.000 Median :0.0000 Median :0.0000 Median :0.00000
Mean :1.847 Mean :0.3915 Mean :0.1958 Mean :0.06349
3rd Qu.:3.000 3rd Qu.:1.0000 3rd Qu.:0.0000 3rd Qu.:0.00000
Max. :3.000 Max. :1.0000 Max. :3.0000 Max. :1.00000
urirr pvft weight agecat
Min. :0.0000 Min. :0.0000 Min. : 709 Min. :1.000
1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:2414 1st Qu.:1.000
Median :0.0000 Median :0.0000 Median :2977 Median :2.000
Mean :0.1481 Mean :0.7937 Mean :2945 Mean :2.238
3rd Qu.:0.0000 3rd Qu.:1.0000 3rd Qu.:3475 3rd Qu.:3.000
Max. :1.0000 Max. :6.0000 Max. :4990 Max. :4.000
wcat anyptl newpvft
Min. :1.000 Min. :0.0000 Min. :0.0000
1st Qu.:2.000 1st Qu.:0.0000 1st Qu.:0.0000
Median :3.000 Median :0.0000 Median :0.0000
Mean :2.857 Mean :0.1587 Mean :0.6931
3rd Qu.:4.000 3rd Qu.:0.0000 3rd Qu.:1.0000
Max. :5.000 Max. :1.0000 Max. :2.0000
# Detailed view of the dataset
glimpse(data.o3)Rows: 189
Columns: 15
$ id <dbl> 85, 86, 87, 88, 89, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 1…
$ lbw <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ age <dbl> 19, 33, 20, 21, 18, 21, 22, 17, 29, 26, 19, 19, 22, 30, 18, 18…
$ lwt <dbl> 182, 155, 105, 108, 107, 124, 118, 103, 123, 113, 95, 150, 95,…
$ race <dbl> 2, 3, 1, 1, 1, 3, 1, 3, 1, 1, 3, 3, 3, 3, 1, 1, 2, 1, 3, 1, 3,…
$ smoke <dbl> 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 1, 0, 1, 0,…
$ ptl <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,…
$ hyper <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,…
$ urirr <dbl> 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0,…
$ pvft <dbl> 0, 3, 1, 2, 0, 0, 1, 1, 1, 0, 0, 1, 0, 2, 0, 0, 0, 3, 0, 1, 2,…
$ weight <dbl> 2523, 2551, 2557, 2594, 2600, 2622, 2637, 2637, 2663, 2665, 27…
$ agecat <dbl+lbl> 1, 4, 2, 2, 1, 2, 2, 1, 3, 3, 1, 1, 2, 4, 1, 1, 1, 3, 2, 3…
$ wcat <dbl+lbl> 5, 5, 1, 2, 2, 3, 2, 1, 3, 2, 1, 4, 1, 2, 1, 1, 1, 2, 2, 2…
$ anyptl <dbl+lbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0…
$ newpvft <dbl+lbl> 0, 2, 1, 2, 0, 0, 1, 1, 1, 0, 0, 1, 0, 2, 0, 0, 0, 2, 0, 1…
4.1.1 Converting Variables to Factor Variables
Since the dataset comes with labeled variables, convert to factor variables using the mutate and across functions from the dplyr package.
# Copy into new dataset and keep original dataset
data.o4 <- data.o3
# Convert to factor variables
data.o4 <- data.o4 %>%
mutate(across(where(is.labelled), ~ as_factor(.x)))
summary(data.o4) id lbw age lwt
Min. : 4.0 Min. :0.0000 Min. :14.00 Min. : 80.0
1st Qu.: 68.0 1st Qu.:0.0000 1st Qu.:19.00 1st Qu.:110.0
Median :123.0 Median :0.0000 Median :23.00 Median :121.0
Mean :121.1 Mean :0.3122 Mean :23.24 Mean :129.8
3rd Qu.:176.0 3rd Qu.:1.0000 3rd Qu.:26.00 3rd Qu.:140.0
Max. :226.0 Max. :1.0000 Max. :45.00 Max. :250.0
race smoke ptl hyper
Min. :1.000 Min. :0.0000 Min. :0.0000 Min. :0.00000
1st Qu.:1.000 1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:0.00000
Median :1.000 Median :0.0000 Median :0.0000 Median :0.00000
Mean :1.847 Mean :0.3915 Mean :0.1958 Mean :0.06349
3rd Qu.:3.000 3rd Qu.:1.0000 3rd Qu.:0.0000 3rd Qu.:0.00000
Max. :3.000 Max. :1.0000 Max. :3.0000 Max. :1.00000
urirr pvft weight agecat wcat
Min. :0.0000 Min. :0.0000 Min. : 709 <20 :51 <105 :37
1st Qu.:0.0000 1st Qu.:0.0000 1st Qu.:2414 20-24:69 106-120:55
Median :0.0000 Median :0.0000 Median :2977 25-29:42 121-130:31
Mean :0.1481 Mean :0.7937 Mean :2945 30+ :27 131-150:30
3rd Qu.:0.0000 3rd Qu.:1.0000 3rd Qu.:3475 151+ :36
Max. :1.0000 Max. :6.0000 Max. :4990
anyptl newpvft
0 :159 0 :100
1+: 30 1 : 47
2+: 42
glimpse(data.o4)Rows: 189
Columns: 15
$ id <dbl> 85, 86, 87, 88, 89, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 1…
$ lbw <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ age <dbl> 19, 33, 20, 21, 18, 21, 22, 17, 29, 26, 19, 19, 22, 30, 18, 18…
$ lwt <dbl> 182, 155, 105, 108, 107, 124, 118, 103, 123, 113, 95, 150, 95,…
$ race <dbl> 2, 3, 1, 1, 1, 3, 1, 3, 1, 1, 3, 3, 3, 3, 1, 1, 2, 1, 3, 1, 3,…
$ smoke <dbl> 0, 0, 1, 1, 1, 0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 1, 0, 1, 0,…
$ ptl <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0,…
$ hyper <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0,…
$ urirr <dbl> 1, 0, 0, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0,…
$ pvft <dbl> 0, 3, 1, 2, 0, 0, 1, 1, 1, 0, 0, 1, 0, 2, 0, 0, 0, 3, 0, 1, 2,…
$ weight <dbl> 2523, 2551, 2557, 2594, 2600, 2622, 2637, 2637, 2663, 2665, 27…
$ agecat <fct> <20, 30+, 20-24, 20-24, <20, 20-24, 20-24, <20, 25-29, 25-29, …
$ wcat <fct> 151+, 151+, <105, 106-120, 106-120, 121-130, 106-120, <105, 12…
$ anyptl <fct> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1+, 0, 0, 0, 0, 0, 0, 0…
$ newpvft <fct> 0, 2+, 1, 2+, 0, 0, 1, 1, 1, 0, 0, 1, 0, 2+, 0, 0, 0, 2+, 0, 1…
4.1.2 Creating a Categorical Outcome Variable
Create a categorical outcome variable using cut function. The variable weight is categorized into a 4-category variable bwt4.
# Create categorical outcome variable
data.o4 <- data.o4 %>%
mutate(bwt4 = cut(weight,
breaks = c(708, 2500, 2999, 3500, max(weight)),
labels = c("<=2500", "2501-3000", "3001-3500", ">3500")
))
# Summarize dataset by the newly created variable `bwt4`
data.o4 %>%
select(bwt4, weight) %>%
group_by(bwt4) %>%
summarize_at(vars(weight), c(min = min, max = max))# A tibble: 4 × 3
bwt4 min max
<fct> <dbl> <dbl>
1 <=2500 709 2495
2 2501-3000 2523 2992
3 3001-3500 3005 3487
4 >3500 3544 4990
4.1.3 Creating an Ordinal Variable
Create an ordinal variable and define the levels of the variable.
# Create and define the levels of ordinal variable
lev <- c(">3500", "3001-3500", "2501-3000", "<=2500")
data.o4 <- data.o4 %>%
mutate(bwt4a = fct_relevel(bwt4, lev)) %>%
mutate(bwt4a = ordered(bwt4a, levels = lev))
# Check the structure of the new variable
str(data.o4$bwt4a) Ord.factor w/ 4 levels ">3500"<"3001-3500"<..: 3 3 3 3 3 3 3 3 3 3 ...
# Check the levels of `bwt4`
levels(data.o4$bwt4)[1] "<=2500" "2501-3000" "3001-3500" ">3500"
# Check the levels of `bwt4a`
levels(data.o4$bwt4a)[1] ">3500" "3001-3500" "2501-3000" "<=2500"
4.2 Baseline Logit Model
The Baseline Logit Model, also known as the Multinomial Logistic Regression Model, is used to handle situations where the outcome variable is nominal with more than two categories. This model extends the binary logistic regression model to multiple categories by modeling the log-odds of each category relative to a baseline category.
This exercise demonstrates how to replicate the baseline logit model (unconstrained) which requires a dataset with an outcome variable without ordering. In this dataset data.o4, a new variable bwt4b without ordering is created from bwt4 by re-leveling the factors according to the specified levels (lev).
# Create new variable and re-leveling according to the specified levels
data.o4 <- data.o4 %>% mutate(bwt4b = fct_relevel(bwt4, lev))
# Fit the model
bl.mod <- vglm(formula = bwt4 ~ smoke, family = multinomial, data = data.o4)
summary(bl.mod)
Call:
vglm(formula = bwt4 ~ smoke, family = multinomial, data = data.o4)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept):1 -0.1881 0.2511 -0.749 0.45392
(Intercept):2 -0.4643 0.2721 -1.707 0.08791 .
(Intercept):3 -0.1881 0.2511 -0.749 0.45392
smoke:1 1.1914 0.4328 2.753 0.00591 **
smoke:2 0.8390 0.4769 1.759 0.07853 .
smoke:3 0.6234 0.4613 1.351 0.17658
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Names of linear predictors: log(mu[,1]/mu[,4]), log(mu[,2]/mu[,4]),
log(mu[,3]/mu[,4])
Residual deviance: 510.9718 on 561 degrees of freedom
Log-likelihood: -255.4859 on 561 degrees of freedom
Number of Fisher scoring iterations: 4
No Hauck-Donner effect found in any of the estimates
Reference group is level 4 of the response
-
smoke:1with an estimate of 1.1914 implies a positive effect of smoking on the probability of being in the first category (<=2500) relative to the reference category (>3500). -
smoke:2andsmoke:3have positive estimates, suggesting an increased likelihood of being in these categories as well due to smoking.
(Intercept):1 (Intercept):2 (Intercept):3 smoke:1 smoke:2
0.8285714 0.6285714 0.8285714 3.2915361 2.3140496
smoke:3
1.8652038
- For
smoke:1, the RRR is 3.2915, meaning smokers are about 3.29 times more likely to be in the lowest birth weight category (<=2500) compared to the reference category (>3500).
4.3 Continuation Ratio Model
The Continuation Ratio Model is another approach to modeling ordinal response data, particularly useful when the response categories have a natural order and when the sequential nature of the response process makes clinical or logical sense. This model is different from other ordinal logistic models in that it models the log-odds of being in a particular category given that the respondent has not been in any previous categories. The model is used to analyze ordinal response data by comparing categories sequentially. In this case, the weight categories are compared in the following manner:
- Compare
bwt >= 3500vs.bwt = 3000-3500 - Compare
bwt >= 3500andbwt = 3000-3500vs.bwt = 2500-3000 - Compare
bwt >= 3500andbwt = 3000-3500andbwt = 2500-3000vs.bwt < 2500
# Check the distribution of weight categories
table(data.o4$bwt4)
<=2500 2501-3000 3001-3500 >3500
59 38 46 46
4.3.1 First Model
# Filter the dataset to include only the weight categories `>3500` and `3001-3500`
data.o4a <- data.o4 %>%
filter(bwt4a == '>3500' | bwt4a == '3001-3500')
# Check the distribution of weight categories
table(data.o4a$bwt4a)
>3500 3001-3500 2501-3000 <=2500
46 46 0 0
# Fit the first model that predicts `bwt4a` based on `smoke`
cr.mod1 <- glm(bwt4a ~ smoke,
family = binomial(link ='logit'),
data = data.o4a)
summary(cr.mod1)
Call:
glm(formula = bwt4a ~ smoke, family = binomial(link = "logit"),
data = data.o4a)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.1881 0.2511 -0.749 0.454
smoke 0.6234 0.4613 1.351 0.177
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 127.54 on 91 degrees of freedom
Residual deviance: 125.68 on 90 degrees of freedom
AIC: 129.68
Number of Fisher Scoring iterations: 4
4.3.2 Second Model
# Filter the dataset to include the weight categories `>3500`, `3001-3500`, and `2501-3000`
data.o4b <- data.o4 %>%
filter(bwt4 == '>3500' | bwt4 == '3001-3500' | bwt4 == '2501-3000')
# Recode `bwt4a` variable into a new variable `bwt4b`
data.o4b <- data.o4b %>%
mutate(bwt4b = ifelse(bwt4a == '>3500', 0, ifelse(bwt4a == '3001-3500', 0, 1)))
# Check the distribution of weight categories
table(data.o4b$bwt4b)
0 1
92 38
# Fit the second model with the recoded `bwt4b` as the dependent variable and `smoke` as the independent variable
cr.mod2 <- glm(bwt4b ~ smoke,
family = binomial(link ='logit'),
data = data.o4b)
summary(cr.mod2)
Call:
glm(formula = bwt4b ~ smoke, family = binomial(link = "logit"),
data = data.o4b)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.0678 0.2471 -4.321 1.56e-05 ***
smoke 0.5082 0.3991 1.273 0.203
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 157.09 on 129 degrees of freedom
Residual deviance: 155.49 on 128 degrees of freedom
AIC: 159.49
Number of Fisher Scoring iterations: 4
4.3.2.1 Third Model
# Recode `bwt4a` variable into `bwt4c` where all categories are grouped into a new variable
data.o4c <- data.o4 %>%
mutate(bwt4c = fct_recode(bwt4a,
gp0 = '>3500',
gp0 = '3001-3500',
gp0 = '3001-3500',
gp0 = '2501-3000'))
# Check the distribution of weight categories
table(data.o4c$bwt4c)
gp0 <=2500
130 59
# Fit the third model to compare the recoded groups with `smoke` as the predictor
cr.mod3 <- glm(bwt4c ~ smoke,
family = binomial(link ='logit'),
data = data.o4c)
summary(cr.mod3)
Call:
glm(formula = bwt4c ~ smoke, family = binomial(link = "logit"),
data = data.o4c)
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.0871 0.2147 -5.062 4.14e-07 ***
smoke 0.7041 0.3196 2.203 0.0276 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 234.67 on 188 degrees of freedom
Residual deviance: 229.80 on 187 degrees of freedom
AIC: 233.8
Number of Fisher Scoring iterations: 4
4.3.3 Conclusion
The three estimates from the different continuation-ratio models are quite similar, with values around 0.6. These estimates indicate that the odds of a birth in the next lower weight category relative to the higher weight categories among women who smoked during pregnancy is about 1.8 times that of women who did not smoke. In summary, smoking during pregnancy consistently increases the odds of having a baby in a lower weight category by approximately 1.8 times across the different comparisons made using the continuation-ratio model.